/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "precomp.hpp"

#ifdef HAVE_VECLIB
#include <vecLib/clapack.h>

typedef __CLPK_integer    integer;
typedef __CLPK_real       real;
#else
#include "clapack.h"
#endif

#undef abs
#undef max
#undef min

namespace cv {

/****************************************************************************************\
*                                 Determinant of the matrix                              *
\****************************************************************************************/

#define det2(m)   (m(0,0)*m(1,1) - m(0,1)*m(1,0))
#define det3(m)   (m(0,0)*(m(1,1)*m(2,2) - m(1,2)*m(2,1)) -  \
                   m(0,1)*(m(1,0)*m(2,2) - m(1,2)*m(2,0)) +  \
                   m(0,2)*(m(1,0)*m(2,1) - m(1,1)*m(2,0)))

double determinant( const Mat& mat ) {
	double result = 0;
	int type = mat.type(), rows = mat.rows;
	size_t step = mat.step;
	const uchar* m = mat.data;

	CV_Assert( mat.rows == mat.cols );

#define Mf(y, x) ((float*)(m + y*step))[x]
#define Md(y, x) ((double*)(m + y*step))[x]

	if ( type == CV_32F ) {
		if ( rows == 2 ) {
			result = det2(Mf);
		} else if ( rows == 3 ) {
			result = det3(Mf);
		} else if ( rows == 1 ) {
			result = Mf(0, 0);
		} else {
			integer i, n = rows, *ipiv, info = 0;
			int bufSize = n * n * sizeof(float) + (n + 1) * sizeof(ipiv[0]), sign = 0;
			AutoBuffer<uchar> buffer(bufSize);

			Mat a(n, n, CV_32F, (uchar*)buffer);
			mat.copyTo(a);

			ipiv = (integer*)cvAlignPtr(a.data + a.step * a.rows, sizeof(integer));
			sgetrf_(&n, &n, (float*)a.data, &n, ipiv, &info);
			assert(info >= 0);

			if ( info == 0 ) {
				result = 1;
				for ( i = 0; i < n; i++ ) {
					result *= ((float*)a.data)[i*(n+1)];
					sign ^= ipiv[i] != i + 1;
				}
				result *= sign ? -1 : 1;
			}
		}
	} else if ( type == CV_64F ) {
		if ( rows == 2 ) {
			result = det2(Md);
		} else if ( rows == 3 ) {
			result = det3(Md);
		} else if ( rows == 1 ) {
			result = Md(0, 0);
		} else {
			integer i, n = rows, *ipiv, info = 0;
			int bufSize = n * n * sizeof(double) + (n + 1) * sizeof(ipiv[0]), sign = 0;
			AutoBuffer<uchar> buffer(bufSize);

			Mat a(n, n, CV_64F, (uchar*)buffer);
			mat.copyTo(a);
			ipiv = (integer*)cvAlignPtr(a.data + a.step * a.rows, sizeof(integer));

			dgetrf_(&n, &n, (double*)a.data, &n, ipiv, &info);
			assert(info >= 0);

			if ( info == 0 ) {
				result = 1;
				for ( i = 0; i < n; i++ ) {
					result *= ((double*)a.data)[i*(n+1)];
					sign ^= ipiv[i] != i + 1;
				}
				result *= sign ? -1 : 1;
			}
		}
	} else {
		CV_Error( CV_StsUnsupportedFormat, "" );
	}

#undef Mf
#undef Md

	return result;
}

/****************************************************************************************\
*                          Inverse (or pseudo-inverse) of a matrix                       *
\****************************************************************************************/

#define Sf( y, x ) ((float*)(srcdata + y*srcstep))[x]
#define Sd( y, x ) ((double*)(srcdata + y*srcstep))[x]
#define Df( y, x ) ((float*)(dstdata + y*dststep))[x]
#define Dd( y, x ) ((double*)(dstdata + y*dststep))[x]

double invert( const Mat& src, Mat& dst, int method ) {
	double result = 0;
	int type = src.type();

	CV_Assert( method == DECOMP_LU || method == DECOMP_CHOLESKY || method == DECOMP_SVD );

	if ( method == DECOMP_SVD ) {
		int n = std::min(src.rows, src.cols);
		SVD svd(src);
		svd.backSubst(Mat(), dst);

		return type == CV_32F ?
			   (((float*)svd.w.data)[0] >= FLT_EPSILON ?
				((float*)svd.w.data)[n-1] / ((float*)svd.w.data)[0] : 0) :
				   (((double*)svd.w.data)[0] >= DBL_EPSILON ?
					((double*)svd.w.data)[n-1] / ((double*)svd.w.data)[0] : 0);
	}

	CV_Assert( src.rows == src.cols && (type == CV_32F || type == CV_64F));
	dst.create( src.rows, src.cols, type );

	if ( method == DECOMP_LU || method == DECOMP_CHOLESKY ) {
		if ( src.rows <= 3 ) {
			uchar* srcdata = src.data;
			uchar* dstdata = dst.data;
			size_t srcstep = src.step;
			size_t dststep = dst.step;

			if ( src.rows == 2 ) {
				if ( type == CV_32FC1 ) {
					double d = det2(Sf);
					if ( d != 0. ) {
						double t0, t1;
						result = d;
						d = 1. / d;
						t0 = Sf(0, 0) * d;
						t1 = Sf(1, 1) * d;
						Df(1, 1) = (float)t0;
						Df(0, 0) = (float)t1;
						t0 = -Sf(0, 1) * d;
						t1 = -Sf(1, 0) * d;
						Df(0, 1) = (float)t0;
						Df(1, 0) = (float)t1;
					}
				} else {
					double d = det2(Sd);
					if ( d != 0. ) {
						double t0, t1;
						result = d;
						d = 1. / d;
						t0 = Sd(0, 0) * d;
						t1 = Sd(1, 1) * d;
						Dd(1, 1) = t0;
						Dd(0, 0) = t1;
						t0 = -Sd(0, 1) * d;
						t1 = -Sd(1, 0) * d;
						Dd(0, 1) = t0;
						Dd(1, 0) = t1;
					}
				}
			} else if ( src.rows == 3 ) {
				if ( type == CV_32FC1 ) {
					double d = det3(Sf);
					if ( d != 0. ) {
						float t[9];
						result = d;
						d = 1. / d;

						t[0] = (float)((Sf(1, 1) * Sf(2, 2) - Sf(1, 2) * Sf(2, 1)) * d);
						t[1] = (float)((Sf(0, 2) * Sf(2, 1) - Sf(0, 1) * Sf(2, 2)) * d);
						t[2] = (float)((Sf(0, 1) * Sf(1, 2) - Sf(0, 2) * Sf(1, 1)) * d);

						t[3] = (float)((Sf(1, 2) * Sf(2, 0) - Sf(1, 0) * Sf(2, 2)) * d);
						t[4] = (float)((Sf(0, 0) * Sf(2, 2) - Sf(0, 2) * Sf(2, 0)) * d);
						t[5] = (float)((Sf(0, 2) * Sf(1, 0) - Sf(0, 0) * Sf(1, 2)) * d);

						t[6] = (float)((Sf(1, 0) * Sf(2, 1) - Sf(1, 1) * Sf(2, 0)) * d);
						t[7] = (float)((Sf(0, 1) * Sf(2, 0) - Sf(0, 0) * Sf(2, 1)) * d);
						t[8] = (float)((Sf(0, 0) * Sf(1, 1) - Sf(0, 1) * Sf(1, 0)) * d);

						Df(0, 0) = t[0]; Df(0, 1) = t[1]; Df(0, 2) = t[2];
						Df(1, 0) = t[3]; Df(1, 1) = t[4]; Df(1, 2) = t[5];
						Df(2, 0) = t[6]; Df(2, 1) = t[7]; Df(2, 2) = t[8];
					}
				} else {
					double d = det3(Sd);
					if ( d != 0. ) {
						double t[9];
						result = d;
						d = 1. / d;

						t[0] = (Sd(1, 1) * Sd(2, 2) - Sd(1, 2) * Sd(2, 1)) * d;
						t[1] = (Sd(0, 2) * Sd(2, 1) - Sd(0, 1) * Sd(2, 2)) * d;
						t[2] = (Sd(0, 1) * Sd(1, 2) - Sd(0, 2) * Sd(1, 1)) * d;

						t[3] = (Sd(1, 2) * Sd(2, 0) - Sd(1, 0) * Sd(2, 2)) * d;
						t[4] = (Sd(0, 0) * Sd(2, 2) - Sd(0, 2) * Sd(2, 0)) * d;
						t[5] = (Sd(0, 2) * Sd(1, 0) - Sd(0, 0) * Sd(1, 2)) * d;

						t[6] = (Sd(1, 0) * Sd(2, 1) - Sd(1, 1) * Sd(2, 0)) * d;
						t[7] = (Sd(0, 1) * Sd(2, 0) - Sd(0, 0) * Sd(2, 1)) * d;
						t[8] = (Sd(0, 0) * Sd(1, 1) - Sd(0, 1) * Sd(1, 0)) * d;

						Dd(0, 0) = t[0]; Dd(0, 1) = t[1]; Dd(0, 2) = t[2];
						Dd(1, 0) = t[3]; Dd(1, 1) = t[4]; Dd(1, 2) = t[5];
						Dd(2, 0) = t[6]; Dd(2, 1) = t[7]; Dd(2, 2) = t[8];
					}
				}
			} else {
				assert( src.rows == 1 );

				if ( type == CV_32FC1 ) {
					double d = Sf(0, 0);
					if ( d != 0. ) {
						result = d;
						Df(0, 0) = (float)(1. / d);
					}
				} else {
					double d = Sd(0, 0);
					if ( d != 0. ) {
						result = d;
						Dd(0, 0) = 1. / d;
					}
				}
			}
			return result;
		}

		src.copyTo(dst);
		integer n = dst.cols, lwork = -1, elem_size = CV_ELEM_SIZE(type),
				lda = (int)(dst.step / elem_size), piv1 = 0, info = 0;

		if ( method == DECOMP_LU ) {
			int buf_size = (int)(n * sizeof(integer));
			AutoBuffer<uchar> buf;
			uchar* buffer;

			if ( type == CV_32F ) {
				real work1 = 0;
				sgetri_(&n, (float*)dst.data, &lda, &piv1, &work1, &lwork, &info);
				lwork = cvRound(work1);
			} else {
				double work1 = 0;
				dgetri_(&n, (double*)dst.data, &lda, &piv1, &work1, &lwork, &info);
				lwork = cvRound(work1);
			}

			buf_size += (int)((lwork + 1) * elem_size);
			buf.allocate(buf_size);
			buffer = (uchar*)buf;

			if ( type == CV_32F ) {
				sgetrf_(&n, &n, (float*)dst.data, &lda, (integer*)buffer, &info);
				if (info == 0)
					sgetri_(&n, (float*)dst.data, &lda, (integer*)buffer,
							(float*)(buffer + n * sizeof(integer)), &lwork, &info);
			} else {
				dgetrf_(&n, &n, (double*)dst.data, &lda, (integer*)buffer, &info);
				if (info == 0)
					dgetri_(&n, (double*)dst.data, &lda, (integer*)buffer,
							(double*)cvAlignPtr(buffer + n * sizeof(integer), elem_size), &lwork, &info);
			}
		} else if ( method == CV_CHOLESKY ) {
			char L[] = {'L', '\0'};
			if ( type == CV_32F ) {
				spotrf_(L, &n, (float*)dst.data, &lda, &info);
				if (info == 0) {
					spotri_(L, &n, (float*)dst.data, &lda, &info);
				}
			} else {
				dpotrf_(L, &n, (double*)dst.data, &lda, &info);
				if (info == 0) {
					dpotri_(L, &n, (double*)dst.data, &lda, &info);
				}
			}
			completeSymm(dst);
		}
		result = info == 0;
	}

	if ( !result ) {
		dst = Scalar(0);
	}

	return result;
}

/****************************************************************************************\
*                              Solving a linear system                                   *
\****************************************************************************************/

bool solve( const Mat& src, const Mat& _src2, Mat& dst, int method ) {
	bool result = true;
	int type = src.type();
	bool is_normal = (method & DECOMP_NORMAL) != 0;

	CV_Assert( type == _src2.type() && (type == CV_32F || type == CV_64F) );

	method &= ~DECOMP_NORMAL;
	CV_Assert( (method != DECOMP_LU && method != DECOMP_CHOLESKY) ||
			   is_normal || src.rows == src.cols );

	// check case of a single equation and small matrix
	if ( (method == DECOMP_LU || method == DECOMP_CHOLESKY) &&
			src.rows <= 3 && src.rows == src.cols && _src2.cols == 1 ) {
		dst.create( src.cols, _src2.cols, src.type() );

#define bf(y) ((float*)(bdata + y*src2step))[0]
#define bd(y) ((double*)(bdata + y*src2step))[0]

		uchar* srcdata = src.data;
		uchar* bdata = _src2.data;
		uchar* dstdata = dst.data;
		size_t srcstep = src.step;
		size_t src2step = _src2.step;
		size_t dststep = dst.step;

		if ( src.rows == 2 ) {
			if ( type == CV_32FC1 ) {
				double d = det2(Sf);
				if ( d != 0. ) {
					float t;
					d = 1. / d;
					t = (float)((bf(0) * Sf(1, 1) - bf(1) * Sf(0, 1)) * d);
					Df(1, 0) = (float)((bf(1) * Sf(0, 0) - bf(0) * Sf(1, 0)) * d);
					Df(0, 0) = t;
				} else {
					result = false;
				}
			} else {
				double d = det2(Sd);
				if ( d != 0. ) {
					double t;
					d = 1. / d;
					t = (bd(0) * Sd(1, 1) - bd(1) * Sd(0, 1)) * d;
					Dd(1, 0) = (bd(1) * Sd(0, 0) - bd(0) * Sd(1, 0)) * d;
					Dd(0, 0) = t;
				} else {
					result = false;
				}
			}
		} else if ( src.rows == 3 ) {
			if ( type == CV_32FC1 ) {
				double d = det3(Sf);
				if ( d != 0. ) {
					float t[3];
					d = 1. / d;

					t[0] = (float)(d *
								   (bf(0) * (Sf(1, 1) * Sf(2, 2) - Sf(1, 2) * Sf(2, 1)) -
									Sf(0, 1) * (bf(1) * Sf(2, 2) - Sf(1, 2) * bf(2)) +
									Sf(0, 2) * (bf(1) * Sf(2, 1) - Sf(1, 1) * bf(2))));

					t[1] = (float)(d *
								   (Sf(0, 0) * (bf(1) * Sf(2, 2) - Sf(1, 2) * bf(2)) -
									bf(0) * (Sf(1, 0) * Sf(2, 2) - Sf(1, 2) * Sf(2, 0)) +
									Sf(0, 2) * (Sf(1, 0) * bf(2) - bf(1) * Sf(2, 0))));

					t[2] = (float)(d *
								   (Sf(0, 0) * (Sf(1, 1) * bf(2) - bf(1) * Sf(2, 1)) -
									Sf(0, 1) * (Sf(1, 0) * bf(2) - bf(1) * Sf(2, 0)) +
									bf(0) * (Sf(1, 0) * Sf(2, 1) - Sf(1, 1) * Sf(2, 0))));

					Df(0, 0) = t[0];
					Df(1, 0) = t[1];
					Df(2, 0) = t[2];
				} else {
					result = false;
				}
			} else {
				double d = det3(Sd);
				if ( d != 0. ) {
					double t[9];

					d = 1. / d;

					t[0] = ((Sd(1, 1) * Sd(2, 2) - Sd(1, 2) * Sd(2, 1)) * bd(0) +
							(Sd(0, 2) * Sd(2, 1) - Sd(0, 1) * Sd(2, 2)) * bd(1) +
							(Sd(0, 1) * Sd(1, 2) - Sd(0, 2) * Sd(1, 1)) * bd(2)) * d;

					t[1] = ((Sd(1, 2) * Sd(2, 0) - Sd(1, 0) * Sd(2, 2)) * bd(0) +
							(Sd(0, 0) * Sd(2, 2) - Sd(0, 2) * Sd(2, 0)) * bd(1) +
							(Sd(0, 2) * Sd(1, 0) - Sd(0, 0) * Sd(1, 2)) * bd(2)) * d;

					t[2] = ((Sd(1, 0) * Sd(2, 1) - Sd(1, 1) * Sd(2, 0)) * bd(0) +
							(Sd(0, 1) * Sd(2, 0) - Sd(0, 0) * Sd(2, 1)) * bd(1) +
							(Sd(0, 0) * Sd(1, 1) - Sd(0, 1) * Sd(1, 0)) * bd(2)) * d;

					Dd(0, 0) = t[0];
					Dd(1, 0) = t[1];
					Dd(2, 0) = t[2];
				} else {
					result = false;
				}
			}
		} else {
			assert( src.rows == 1 );

			if ( type == CV_32FC1 ) {
				double d = Sf(0, 0);
				if ( d != 0. ) {
					Df(0, 0) = (float)(bf(0) / d);
				} else {
					result = false;
				}
			} else {
				double d = Sd(0, 0);
				if ( d != 0. ) {
					Dd(0, 0) = (bd(0) / d);
				} else {
					result = false;
				}
			}
		}
		return result;
	}

	double rcond = -1, s1 = 0, work1 = 0, *work = 0, *s = 0;
	float frcond = -1, fs1 = 0, fwork1 = 0, *fwork = 0, *fs = 0;
	integer m = src.rows, m_ = m, n = src.cols, mn = std::max(m, n),
			nm = std::min(m, n), nb = _src2.cols, lwork = -1, liwork = 0, iwork1 = 0,
			lda = m, ldx = mn, info = 0, rank = 0, *iwork = 0;
	int elem_size = CV_ELEM_SIZE(type);
	bool copy_rhs = false;
	int buf_size = 0;
	AutoBuffer<uchar> buffer;
	uchar* ptr;
	char N[] = {'N', '\0'}, L[] = {'L', '\0'};

	Mat src2 = _src2;
	dst.create( src.cols, src2.cols, src.type() );

	if ( m <= n ) {
		is_normal = false;
	} else if ( is_normal ) {
		m_ = n;
	}

	buf_size += (is_normal ? n * n : m * n) * elem_size;

	if ( m_ != n || nb > 1 || !dst.isContinuous() ) {
		copy_rhs = true;
		if ( is_normal ) {
			buf_size += n * nb * elem_size;
		} else {
			buf_size += mn * nb * elem_size;
		}
	}

	if ( method == DECOMP_SVD || method == DECOMP_EIG ) {
		integer nlvl = cvRound(std::log(std::max(std::min(m_, n) / 25., 1.)) / CV_LOG2) + 1;
		liwork = std::min(m_, n) * (3 * std::max(nlvl, (integer)0) + 11);

		if ( type == CV_32F )
			sgelsd_(&m_, &n, &nb, (float*)src.data, &lda, (float*)dst.data, &ldx,
					&fs1, &frcond, &rank, &fwork1, &lwork, &iwork1, &info);
		else
			dgelsd_(&m_, &n, &nb, (double*)src.data, &lda, (double*)dst.data, &ldx,
					&s1, &rcond, &rank, &work1, &lwork, &iwork1, &info );
		buf_size += nm * elem_size + (liwork + 1) * sizeof(integer);
	} else if ( method == DECOMP_QR ) {
		if ( type == CV_32F )
			sgels_(N, &m_, &n, &nb, (float*)src.data, &lda,
				   (float*)dst.data, &ldx, &fwork1, &lwork, &info );
		else
			dgels_(N, &m_, &n, &nb, (double*)src.data, &lda,
				   (double*)dst.data, &ldx, &work1, &lwork, &info );
	} else if ( method == DECOMP_LU ) {
		buf_size += (n + 1) * sizeof(integer);
	} else if ( method == DECOMP_CHOLESKY )
		;
	else {
		CV_Error( CV_StsBadArg, "Unknown method" );
	}
	assert(info == 0);

	lwork = cvRound(type == CV_32F ? (double)fwork1 : work1);
	buf_size += lwork * elem_size;
	buffer.allocate(buf_size);
	ptr = (uchar*)buffer;

	Mat at(n, m_, type, ptr);
	ptr += n * m_ * elem_size;

	if ( method == DECOMP_CHOLESKY || method == DECOMP_EIG ) {
		src.copyTo(at);
	} else if ( !is_normal ) {
		transpose(src, at);
	} else {
		mulTransposed(src, at, true);
	}

	Mat xt;
	if ( !is_normal ) {
		if ( copy_rhs ) {
			Mat temp(nb, mn, type, ptr);
			ptr += nb * mn * elem_size;
			Mat bt = temp.colRange(0, m);
			xt = temp.colRange(0, n);
			transpose(src2, bt);
		} else {
			src2.copyTo(dst);
			xt = Mat(1, n, type, dst.data);
		}
	} else {
		if ( copy_rhs ) {
			xt = Mat(nb, n, type, ptr);
			ptr += nb * n * elem_size;
		} else {
			xt = Mat(1, n, type, dst.data);
		}
		// (a'*b)' = b'*a
		gemm( src2, src, 1, Mat(), 0, xt, GEMM_1_T );
	}

	lda = (int)(at.step ? at.step / elem_size : at.cols);
	ldx = (int)(xt.step ? xt.step / elem_size : (!is_normal && copy_rhs ? mn : n));

	if ( method == DECOMP_SVD || method == DECOMP_EIG ) {
		if ( type == CV_32F ) {
			fs = (float*)ptr;
			ptr += nm * elem_size;
			fwork = (float*)ptr;
			ptr += lwork * elem_size;
			iwork = (integer*)cvAlignPtr(ptr, sizeof(integer));

			sgelsd_(&m_, &n, &nb, (float*)at.data, &lda, (float*)xt.data, &ldx,
					fs, &frcond, &rank, fwork, &lwork, iwork, &info);
		} else {
			s = (double*)ptr;
			ptr += nm * elem_size;
			work = (double*)ptr;
			ptr += lwork * elem_size;
			iwork = (integer*)cvAlignPtr(ptr, sizeof(integer));

			dgelsd_(&m_, &n, &nb, (double*)at.data, &lda, (double*)xt.data, &ldx,
					s, &rcond, &rank, work, &lwork, iwork, &info);
		}
	} else if ( method == CV_QR ) {
		if ( type == CV_32F ) {
			fwork = (float*)ptr;
			sgels_(N, &m_, &n, &nb, (float*)at.data, &lda,
				   (float*)xt.data, &ldx, fwork, &lwork, &info);
		} else {
			work = (double*)ptr;
			dgels_(N, &m_, &n, &nb, (double*)at.data, &lda,
				   (double*)xt.data, &ldx, work, &lwork, &info);
		}
	} else if ( method == CV_CHOLESKY || (method == CV_LU && is_normal) ) {
		if ( type == CV_32F ) {
			spotrf_(L, &n, (float*)at.data, &lda, &info);
			if (info == 0) {
				spotrs_(L, &n, &nb, (float*)at.data, &lda, (float*)xt.data, &ldx, &info);
			}
		} else {
			dpotrf_(L, &n, (double*)at.data, &lda, &info);
			if (info == 0) {
				dpotrs_(L, &n, &nb, (double*)at.data, &lda, (double*)xt.data, &ldx, &info);
			}
		}
	} else if ( method == CV_LU ) {
		iwork = (integer*)cvAlignPtr(ptr, sizeof(integer));
		if ( type == CV_32F ) {
			sgesv_(&n, &nb, (float*)at.data, &lda, iwork, (float*)xt.data, &ldx, &info );
		} else {
			dgesv_(&n, &nb, (double*)at.data, &lda, iwork, (double*)xt.data, &ldx, &info );
		}
	} else {
		assert(0);
	}
	result = info == 0;

	if ( !result ) {
		dst = Scalar(0);
	} else if ( xt.data != dst.data ) {
		transpose( xt, dst );
	}

	return result;
}


/////////////////// finding eigenvalues and eigenvectors of a symmetric matrix ///////////////

template<typename Real> static inline Real hypot(Real a, Real b) {
	a = std::abs(a);
	b = std::abs(b);
	Real f;
	if ( a > b ) {
		f = b / a;
		return a * std::sqrt(1 + f * f);
	}
	if ( b == 0 ) {
		return 0;
	}
	f = a / b;
	return b * std::sqrt(1 + f * f);
}


template<typename Real> bool jacobi(const Mat& _S0, Mat& _e, Mat& matE, bool computeEvects, Real eps) {
	int n = _S0.cols, i, j, k, m;

	if ( computeEvects ) {
		matE = Mat::eye(n, n, _S0.type());
	}

	int iters, maxIters = n * n * 30;

	AutoBuffer<uchar> buf(n * 2 * sizeof(int) + (n * n + n * 2 + 1)*sizeof(Real));
	Real* S = alignPtr((Real*)(uchar*)buf, sizeof(Real));
	Real* maxSR = S + n * n;
	Real* maxSC = maxSR + n;
	int* indR = (int*)(maxSC + n);
	int* indC = indR + n;

	Mat matS(_S0.size(), _S0.type(), S);
	_S0.copyTo(matS);

	Real mv;
	Real* E = (Real*)matE.data;
	Real* e = (Real*)_e.data;
	int Sstep = matS.step / sizeof(Real);
	int estep = _e.rows == 1 ? 1 : _e.step / sizeof(Real);
	int Estep = matE.step / sizeof(Real);

	for ( k = 0; k < n; k++ ) {
		e[k* estep] = S[(Sstep + 1)*k];
		if ( k < n - 1 ) {
			for ( m = k + 1, mv = std::abs(S[Sstep*k + m]), i = k + 2; i < n; i++ ) {
				Real v = std::abs(S[Sstep*k+i]);
				if ( mv < v ) {
					mv = v, m = i;
				}
			}
			maxSR[k] = mv;
			indR[k] = m;
		}
		if ( k > 0 ) {
			for ( m = 0, mv = std::abs(S[k]), i = 1; i < k; i++ ) {
				Real v = std::abs(S[Sstep*i+k]);
				if ( mv < v ) {
					mv = v, m = i;
				}
			}
			maxSC[k] = mv;
			indC[k] = m;
		}
	}

	for ( iters = 0; iters < maxIters; iters++ ) {
		// find index (k,l) of pivot p
		for ( k = 0, mv = maxSR[0], i = 1; i < n - 1; i++ ) {
			Real v = maxSR[i];
			if ( mv < v ) {
				mv = v, k = i;
			}
		}
		int l = indR[k];
		for ( i = 1; i < n; i++ ) {
			Real v = maxSC[i];
			if ( mv < v ) {
				mv = v, k = indC[i], l = i;
			}
		}

		Real p = S[Sstep*k + l];
		if ( std::abs(p) <= eps ) {
			break;
		}
		Real y = Real((e[estep*l] - e[estep*k]) * 0.5);
		Real t = std::abs(y) + hypot(p, y);
		Real s = hypot(p, t);
		Real c = t / s;
		s = p / s; t = (p / t) * p;
		if ( y < 0 ) {
			s = -s, t = -t;
		}
		S[Sstep* k + l] = 0;

		e[estep* k] -= t;
		e[estep* l] += t;

		Real a0, b0;

#undef rotate
#define rotate(v0, v1) a0 = v0, b0 = v1, v0 = a0*c - b0*s, v1 = a0*s + b0*c

		// rotate rows and columns k and l
		for ( i = 0; i < k; i++ ) {
			rotate(S[Sstep*i+k], S[Sstep*i+l]);
		}
		for ( i = k + 1; i < l; i++ ) {
			rotate(S[Sstep*k+i], S[Sstep*i+l]);
		}
		for ( i = l + 1; i < n; i++ ) {
			rotate(S[Sstep*k+i], S[Sstep*l+i]);
		}

		// rotate eigenvectors
		if ( computeEvects )
			for ( i = 0; i < n; i++ ) {
				rotate(E[Estep*k+i], E[Estep*l+i]);
			}

#undef rotate

		for ( j = 0; j < 2; j++ ) {
			int idx = j == 0 ? k : l;
			if ( idx < n - 1 ) {
				for ( m = idx + 1, mv = std::abs(S[Sstep*idx + m]), i = idx + 2; i < n; i++ ) {
					Real v = std::abs(S[Sstep*idx+i]);
					if ( mv < v ) {
						mv = v, m = i;
					}
				}
				maxSR[idx] = mv;
				indR[idx] = m;
			}
			if ( idx > 0 ) {
				for ( m = 0, mv = std::abs(S[idx]), i = 1; i < idx; i++ ) {
					Real v = std::abs(S[Sstep*i+idx]);
					if ( mv < v ) {
						mv = v, m = i;
					}
				}
				maxSC[idx] = mv;
				indC[idx] = m;
			}
		}
	}

	// sort eigenvalues & eigenvectors
	for ( k = 0; k < n - 1; k++ ) {
		m = k;
		for ( i = k + 1; i < n; i++ ) {
			if ( e[estep*m] < e[estep*i] ) {
				m = i;
			}
		}
		if ( k != m ) {
			std::swap(e[estep*m], e[estep*k]);
			if ( computeEvects )
				for ( i = 0; i < n; i++ ) {
					std::swap(E[Estep*m + i], E[Estep*k + i]);
				}
		}
	}

	return true;
}


static bool eigen( const Mat& src, Mat& evals, Mat& evects, bool computeEvects,
				   int lowindex, int highindex ) {
	int type = src.type();
	integer n = src.rows;

	// If a range is selected both limits are needed.
	CV_Assert( ( lowindex >= 0 && highindex >= 0 ) ||
			   ( lowindex < 0 && highindex < 0 ) );

	// lapack sorts from lowest to highest so we flip
	integer il = n - highindex;
	integer iu = n - lowindex;

	CV_Assert( src.rows == src.cols );
	CV_Assert (type == CV_32F || type == CV_64F);

	// allow for 1xn eigenvalue matrix too
	if ( !(evals.rows == 1 && evals.cols == n && evals.type() == type) ) {
		evals.create(n, 1, type);
	}

	if ( n <= 20 ) {
		if ( type == CV_32F ) {
			return jacobi<float>(src, evals, evects, computeEvects, FLT_EPSILON);
		} else {
			return jacobi<double>(src, evals, evects, computeEvects, DBL_EPSILON);
		}
	}

	bool result;
	integer m = 0, lda, ldv = n, lwork = -1, iwork1 = 0, liwork = -1, idummy = 0, info = 0;
	integer* isupport, *iwork;
	char job[] = { computeEvects ? 'V' : 'N', '\0' };
	char range[2] = "I";
	range[0] = (il < n + 1) ? 'I' : 'A';

	char L[] = {'L', '\0'};
	uchar* work;

	AutoBuffer<uchar> buf;

	int elem_size = (int)src.elemSize();
	lda = (int)(src.step / elem_size);

	if ( computeEvects ) {
		evects.create(n, n, type);
		ldv = (int)(evects.step / elem_size);
	}

	bool copy_evals = !evals.isContinuous();

	if ( type == CV_32FC1 ) {
		float work1 = 0, dummy = 0, abstol = 0, *s;

		ssyevr_(job, range, L, &n, (float*)src.data, &lda, &dummy, &dummy, &il, &iu,
				&abstol, &m, (float*)evals.data, (float*)evects.data, &ldv,
				&idummy, &work1, &lwork, &iwork1, &liwork, &info );
		assert( info == 0 );

		lwork = cvRound(work1);
		liwork = iwork1;
		buf.allocate((lwork + n * n + (copy_evals ? n : 0))*elem_size +
					 (liwork + 2 * n + 1)*sizeof(integer));
		Mat a(n, n, type, (uchar*)buf);
		src.copyTo(a);
		lda = a.step1();
		work = a.data + n * n * elem_size;
		if ( copy_evals ) {
			s = (float*)(work + lwork * elem_size);
		} else {
			s = (float*)evals.data;
		}

		iwork = (integer*)cvAlignPtr(work + (lwork + (copy_evals ? n : 0)) * elem_size, sizeof(integer));
		isupport = iwork + liwork;

		ssyevr_(job, range, L, &n, (float*)a.data, &lda, &dummy, &dummy,
				&il, &iu, &abstol, &m, s, (float*)evects.data,
				&ldv, isupport, (float*)work, &lwork, iwork, &liwork, &info );
		result = info == 0;
	} else {
		double work1 = 0, dummy = 0, abstol = 0, *s;

		dsyevr_(job, range, L, &n, (double*)src.data, &lda, &dummy, &dummy, &il, &iu,
				&abstol, &m, (double*)evals.data, (double*)evects.data, &ldv,
				&idummy, &work1, &lwork, &iwork1, &liwork, &info );
		assert( info == 0 );

		lwork = cvRound(work1);
		liwork = iwork1;
		buf.allocate((lwork + n * n + (copy_evals ? n : 0))*elem_size +
					 (liwork + 2 * n + 1)*sizeof(integer));
		Mat a(n, n, type, (uchar*)buf);
		src.copyTo(a);
		lda = a.step1();
		work = a.data + n * n * elem_size;

		if ( copy_evals ) {
			s = (double*)(work + lwork * elem_size);
		} else {
			s = (double*)evals.data;
		}

		iwork = (integer*)cvAlignPtr(work + (lwork + (copy_evals ? n : 0)) * elem_size, sizeof(integer));
		isupport = iwork + liwork;

		dsyevr_(job, range, L, &n, (double*)a.data, &lda, &dummy, &dummy,
				&il, &iu, &abstol, &m, s, (double*)evects.data,
				&ldv, isupport, (double*)work, &lwork, iwork, &liwork, &info );
		result = info == 0;
	}

	if ( copy_evals ) {
		Mat(evals.rows, evals.cols, type, work + lwork * elem_size).copyTo(evals);
	}

	if ( il < n + 1 && n > 20 ) {
		int nVV = iu - il + 1;
		if ( computeEvects ) {
			Mat flipme = evects.rowRange(0, nVV);
			flip(flipme, flipme, 0);
			flipme = evals.rowRange(0, nVV);
			flip(flipme, flipme, 0);
		}
	} else {
		flip(evals, evals, evals.rows > 1 ? 0 : 1);
		if ( computeEvects ) {
			flip(evects, evects, 0);
		}
	}

	return result;
}

bool eigen( const Mat& src, Mat& evals, int lowindex, int highindex ) {
	Mat evects;
	return eigen(src, evals, evects, false, lowindex, highindex);
}

bool eigen( const Mat& src, Mat& evals, Mat& evects, int lowindex,
			int highindex ) {
	return eigen(src, evals, evects, true, lowindex, highindex);
}



/* y[0:m,0:n] += diag(a[0:1,0:m]) * x[0:m,0:n] */
template<typename T1, typename T2, typename T3> static void
MatrAXPY( int m, int n, const T1* x, int dx,
		  const T2* a, int inca, T3* y, int dy ) {
	int i, j;
	for ( i = 0; i < m; i++, x += dx, y += dy ) {
		T2 s = a[i*inca];
		for ( j = 0; j <= n - 4; j += 4 ) {
			T3 t0 = (T3)(y[j]   + s * x[j]);
			T3 t1 = (T3)(y[j+1] + s * x[j+1]);
			y[j]   = t0;
			y[j+1] = t1;
			t0 = (T3)(y[j+2] + s * x[j+2]);
			t1 = (T3)(y[j+3] + s * x[j+3]);
			y[j+2] = t0;
			y[j+3] = t1;
		}

		for ( ; j < n; j++ ) {
			y[j] = (T3)(y[j] + s * x[j]);
		}
	}
}

template<typename T> static void
SVBkSb( int m, int n, const T* w, int incw,
		const T* u, int ldu, int uT,
		const T* v, int ldv, int vT,
		const T* b, int ldb, int nb,
		T* x, int ldx, double* buffer, T eps ) {
	double threshold = 0;
	int udelta0 = uT ? ldu : 1, udelta1 = uT ? 1 : ldu;
	int vdelta0 = vT ? ldv : 1, vdelta1 = vT ? 1 : ldv;
	int i, j, nm = std::min(m, n);

	if ( !b ) {
		nb = m;
	}

	for ( i = 0; i < n; i++ )
		for ( j = 0; j < nb; j++ ) {
			x[i* ldx + j] = 0;
		}

	for ( i = 0; i < nm; i++ ) {
		threshold += w[i*incw];
	}
	threshold *= eps;

	// v * inv(w) * uT * b
	for ( i = 0; i < nm; i++, u += udelta0, v += vdelta0 ) {
		double wi = w[i*incw];
		if ( wi <= threshold ) {
			continue;
		}
		wi = 1 / wi;

		if ( nb == 1 ) {
			double s = 0;
			if ( b )
				for ( j = 0; j < m; j++ ) {
					s += u[j*udelta1] * b[j*ldb];
				}
			else {
				s = u[0];
			}
			s *= wi;

			for ( j = 0; j < n; j++ ) {
				x[j* ldx] = (T)(x[j*ldx] + s * v[j*vdelta1]);
			}
		} else {
			if ( b ) {
				for ( j = 0; j < nb; j++ ) {
					buffer[j] = 0;
				}
				MatrAXPY( m, nb, b, ldb, u, udelta1, buffer, 0 );
				for ( j = 0; j < nb; j++ ) {
					buffer[j] *= wi;
				}
			} else {
				for ( j = 0; j < nb; j++ ) {
					buffer[j] = u[j*udelta1] * wi;
				}
			}
			MatrAXPY( n, nb, buffer, 0, v, vdelta1, x, ldx );
		}
	}
}


SVD& SVD::operator ()(const Mat& a, int flags) {
	integer m = a.rows, n = a.cols, mn = std::max(m, n), nm = std::min(m, n);
	int type = a.type(), elem_size = (int)a.elemSize();

	if ( flags & NO_UV ) {
		u.release();
		vt.release();
	} else {
		u.create( (int)m, (int)((flags & FULL_UV) ? m : nm), type );
		vt.create( (int)((flags & FULL_UV) ? n : nm), n, type );
	}

	w.create(nm, 1, type);

	Mat _a = a;
	int a_ofs = 0, work_ofs = 0, iwork_ofs = 0, buf_size = 0;
	bool temp_a = false;
	double u1 = 0, v1 = 0, work1 = 0;
	float uf1 = 0, vf1 = 0, workf1 = 0;
	integer lda, ldu, ldv, lwork = -1, iwork1 = 0, info = 0, *iwork = 0;
	char mode[] = {u.data || vt.data ? 'S' : 'N', '\0'};

	if ( m != n && !(flags & NO_UV) && (flags & FULL_UV) ) {
		mode[0] = 'A';
	}

	if ( !(flags & MODIFY_A) ) {
		if ( mode[0] == 'N' || mode[0] == 'A' ) {
			temp_a = true;
		} else if ( ((vt.data && a.size() == vt.size()) || (u.data && a.size() == u.size())) &&
					mode[0] == 'S' ) {
			mode[0] = 'O';
		}
	}

	lda = a.cols;
	ldv = ldu = mn;

	if ( type == CV_32F ) {
		sgesdd_(mode, &n, &m, (float*)a.data, &lda, (float*)w.data,
				&vf1, &ldv, &uf1, &ldu, &workf1, &lwork, &iwork1, &info );
		lwork = cvRound(workf1);
	} else {
		dgesdd_(mode, &n, &m, (double*)a.data, &lda, (double*)w.data,
				&v1, &ldv, &u1, &ldu, &work1, &lwork, &iwork1, &info );
		lwork = cvRound(work1);
	}

	assert(info == 0);
	if ( temp_a ) {
		a_ofs = buf_size;
		buf_size += n * m * elem_size;
	}
	work_ofs = buf_size;
	buf_size += lwork * elem_size;
	buf_size = cvAlign(buf_size, sizeof(iwork[0]));
	iwork_ofs = buf_size;
	buf_size += 8 * nm * sizeof(integer);

	AutoBuffer<uchar> buf(buf_size);
	uchar* buffer = (uchar*)buf;

	if ( temp_a ) {
		_a = Mat(a.rows, a.cols, type, buffer );
		a.copyTo(_a);
	}

	if ( !(flags & MODIFY_A) && !temp_a ) {
		if ( vt.data && a.size() == vt.size() ) {
			a.copyTo(vt);
			_a = vt;
		} else if ( u.data && a.size() == u.size() ) {
			a.copyTo(u);
			_a = u;
		}
	}

	if ( mode[0] != 'N' ) {
		ldv = (int)(vt.step ? vt.step / elem_size : vt.cols);
		ldu = (int)(u.step ? u.step / elem_size : u.cols);
	}

	lda = (int)(_a.step ? _a.step / elem_size : _a.cols);
	if ( type == CV_32F ) {
		sgesdd_(mode, &n, &m, (float*)_a.data, &lda, (float*)w.data,
				(float*)vt.data, &ldv, (float*)u.data, &ldu,
				(float*)(buffer + work_ofs), &lwork, (integer*)(buffer + iwork_ofs), &info );
	} else {
		dgesdd_(mode, &n, &m, (double*)_a.data, &lda, (double*)w.data,
				(double*)vt.data, &ldv, (double*)u.data, &ldu,
				(double*)(buffer + work_ofs), &lwork, (integer*)(buffer + iwork_ofs), &info );
	}
	CV_Assert(info >= 0);
	if (info != 0) {
		u = Scalar(0.);
		vt = Scalar(0.);
		w = Scalar(0.);
	}
	return *this;
}


void SVD::backSubst( const Mat& rhs, Mat& dst ) const {
	int type = w.type(), esz = (int)w.elemSize();
	int m = u.rows, n = vt.cols, nb = rhs.data ? rhs.cols : m;
	AutoBuffer<double> buffer(nb);
	CV_Assert( u.data && vt.data && w.data );

	if ( rhs.data ) {
		CV_Assert( rhs.type() == type && rhs.rows == m );
	}

	dst.create( n, nb, type );
	if ( type == CV_32F )
		SVBkSb(m, n, (float*)w.data, 1, (float*)u.data, (int)(u.step / esz), false,
			   (float*)vt.data, (int)(vt.step / esz), true, (float*)rhs.data, (int)(rhs.step / esz),
			   nb, (float*)dst.data, (int)(dst.step / esz), buffer, 10 * FLT_EPSILON );
	else if ( type == CV_64F )
		SVBkSb(m, n, (double*)w.data, 1, (double*)u.data, (int)(u.step / esz), false,
			   (double*)vt.data, (int)(vt.step / esz), true, (double*)rhs.data, (int)(rhs.step / esz),
			   nb, (double*)dst.data, (int)(dst.step / esz), buffer, 2 * DBL_EPSILON );
	else {
		CV_Error( CV_StsUnsupportedFormat, "" );
	}
}

}


CV_IMPL double
cvDet( const CvArr* arr ) {
	if ( CV_IS_MAT(arr) && ((CvMat*)arr)->rows <= 3 ) {
		CvMat* mat = (CvMat*)arr;
		int type = CV_MAT_TYPE(mat->type);
		int rows = mat->rows;
		uchar* m = mat->data.ptr;
		int step = mat->step;
		CV_Assert( rows == mat->cols );

#define Mf(y, x) ((float*)(m + y*step))[x]
#define Md(y, x) ((double*)(m + y*step))[x]

		if ( type == CV_32F ) {
			if ( rows == 2 ) {
				return det2(Mf);
			}
			if ( rows == 3 ) {
				return det3(Mf);
			}
		} else if ( type == CV_64F ) {
			if ( rows == 2 ) {
				return det2(Md);
			}
			if ( rows == 3 ) {
				return det3(Md);
			}
		}
		return cv::determinant(cv::Mat(mat));
	}
	return cv::determinant(cv::cvarrToMat(arr));
}


CV_IMPL double
cvInvert( const CvArr* srcarr, CvArr* dstarr, int method ) {
	cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);

	CV_Assert( src.type() == dst.type() && src.rows == dst.cols && src.cols == dst.rows );
	return cv::invert( src, dst, method == CV_CHOLESKY ? cv::DECOMP_CHOLESKY :
					   method == CV_SVD || method == CV_SVD_SYM ? cv::DECOMP_SVD : cv::DECOMP_LU );
}


CV_IMPL int
cvSolve( const CvArr* Aarr, const CvArr* barr, CvArr* xarr, int method ) {
	cv::Mat A = cv::cvarrToMat(Aarr), b = cv::cvarrToMat(barr), x = cv::cvarrToMat(xarr);

	CV_Assert( A.type() == x.type() && A.cols == x.rows && x.cols == b.cols );
	return cv::solve( A, b, x, method == CV_CHOLESKY ? cv::DECOMP_CHOLESKY :
					  method == CV_SVD || method == CV_SVD_SYM ? cv::DECOMP_SVD :
					  A.rows > A.cols ? cv::DECOMP_QR : cv::DECOMP_LU );
}


CV_IMPL void
cvEigenVV( CvArr* srcarr, CvArr* evectsarr, CvArr* evalsarr, double,
		   int lowindex, int highindex) {
	cv::Mat src = cv::cvarrToMat(srcarr), evals = cv::cvarrToMat(evalsarr);
	if ( evectsarr ) {
		cv::Mat evects = cv::cvarrToMat(evectsarr);
		eigen(src, evals, evects, lowindex, highindex);
	} else {
		eigen(src, evals, lowindex, highindex);
	}
}


CV_IMPL void
cvSVD( CvArr* aarr, CvArr* warr, CvArr* uarr, CvArr* varr, int flags ) {
	cv::Mat a = cv::cvarrToMat(aarr), w = cv::cvarrToMat(warr), u, v;
	int m = a.rows, n = a.cols, type = a.type(), mn = std::max(m, n), nm = std::min(m, n);

	CV_Assert( w.type() == type &&
			   (w.size() == cv::Size(nm, 1) || w.size() == cv::Size(1, nm) ||
				w.size() == cv::Size(nm, nm) || w.size() == cv::Size(n, m)) );

	cv::SVD svd;

	if ( w.size() == cv::Size(nm, 1) ) {
		svd.w = cv::Mat(nm, 1, type, w.data );
	} else if ( w.isContinuous() ) {
		svd.w = w;
	}

	if ( uarr ) {
		u = cv::cvarrToMat(uarr);
		CV_Assert( u.type() == type );
		svd.u = u;
	}

	if ( varr ) {
		v = cv::cvarrToMat(varr);
		CV_Assert( v.type() == type );
		svd.vt = v;
	}

	svd(a, ((flags & CV_SVD_MODIFY_A) ? cv::SVD::MODIFY_A : 0) |
		((!svd.u.data && !svd.vt.data) ? cv::SVD::NO_UV : 0) |
		((m != n && (svd.u.size() == cv::Size(mn, mn) ||
					 svd.vt.size() == cv::Size(mn, mn))) ? cv::SVD::FULL_UV : 0));

	if ( u.data ) {
		if ( flags & CV_SVD_U_T ) {
			cv::transpose( svd.u, u );
		} else if ( u.data != svd.u.data ) {
			CV_Assert( u.size() == svd.u.size() );
			svd.u.copyTo(u);
		}
	}

	if ( v.data ) {
		if ( !(flags & CV_SVD_V_T) ) {
			cv::transpose( svd.vt, v );
		} else if ( v.data != svd.vt.data ) {
			CV_Assert( v.size() == svd.vt.size() );
			svd.vt.copyTo(v);
		}
	}

	if ( w.data != svd.w.data ) {
		if ( w.size() == svd.w.size() ) {
			svd.w.copyTo(w);
		} else {
			w = cv::Scalar(0);
			cv::Mat wd = w.diag();
			svd.w.copyTo(wd);
		}
	}
}


CV_IMPL void
cvSVBkSb( const CvArr* warr, const CvArr* uarr,
		  const CvArr* varr, const CvArr* rhsarr,
		  CvArr* dstarr, int flags ) {
	cv::Mat w = cv::cvarrToMat(warr), u = cv::cvarrToMat(uarr),
			v = cv::cvarrToMat(varr), rhs, dst = cv::cvarrToMat(dstarr);
	int type = w.type();
	bool uT = (flags & CV_SVD_U_T) != 0, vT = (flags & CV_SVD_V_T) != 0;
	int m = !uT ? u.rows : u.cols;
	int n = vT ? v.cols : v.rows;
	int nm = std::min(n, m), nb;
	int esz = (int)w.elemSize();
	int incw = w.size() == cv::Size(nm, 1) ? 1 : (int)(w.step / esz) + (w.cols > 1 && w.rows > 1);

	CV_Assert( type == u.type() && type == v.type() &&
			   type == dst.type() && dst.rows == n &&
			   (!uT ? u.cols : u.rows) >= nm && (vT ? v.rows : v.cols) >= nm &&
			   (w.size() == cv::Size(nm, 1) || w.size() == cv::Size(1, nm) ||
				w.size() == cv::Size(nm, nm) || w.size() == cv::Size(n, m)));

	if ( rhsarr ) {
		rhs = cv::cvarrToMat(rhsarr);
		nb = rhs.cols;
		CV_Assert( type == rhs.type() );
	} else {
		nb = m;
	}

	CV_Assert( dst.cols == nb );
	cv::AutoBuffer<double> buffer(nb);

	if ( type == CV_32F )
		cv::SVBkSb(m, n, (float*)w.data, incw, (float*)u.data, (int)(u.step / esz), uT,
				   (float*)v.data, (int)(v.step / esz), vT, (float*)rhs.data, (int)(rhs.step / esz),
				   nb, (float*)dst.data, (int)(dst.step / esz), buffer, 2 * FLT_EPSILON );
	else
		cv::SVBkSb(m, n, (double*)w.data, incw, (double*)u.data, (int)(u.step / esz), uT,
				   (double*)v.data, (int)(v.step / esz), vT, (double*)rhs.data, (int)(rhs.step / esz),
				   nb, (double*)dst.data, (int)(dst.step / esz), buffer, 2 * DBL_EPSILON );
}
